library(foreign)

genperm <- function(){
	return(c(as.vector(replicate(63,sample(c(0,0,1),3))),sample(c(0,0,1),2)))
}

#RI Function for ATE estimation
est.ate <- function(treat, outcome){
	return(mean(outcome[treat==1]) - mean(outcome[treat==0]))
}

ri <- function(outcome, treat, perms){
	ate <- est.ate(treat, outcome)
	ate.dist.under.sharp.null <- apply(perms, 2, est.ate, outcome)
	p.value <- mean(ate <= ate.dist.under.sharp.null)
	plot(density(ate.dist.under.sharp.null))
	abline(v = ate)
	return(list(ate = ate, p.value = p.value))
}

#Read in overall data
data <- read.dta("kalla-broockman-donor-access-2013-data.dta", convert.underscore = TRUE)
data <- data[order(data$block),] #put data in order of blocks to make RI easier
data <- data[c(1:135,138:191,136,137),] #put block 46 at end since it only has two in it
head(data)


#Results
table(data$treat.donor, data$staffrank)

#Table 1 point estimates
round(table(data$treat.donor, data$staffrank)/rowSums(table(data$treat.donor, data$staffrank)), digits = 3)*100

#Cumulative probs
cumulative.probs <- matrix(NA, nrow = 2, ncol = 6)
for(i in 6:1){
	for(j in 1:2){
		cumulative.probs[j,7-i] <- sum(table(data$treat.donor, data$staffrank)[j,i:6])
	}
}
cumulative.probs[1,] <- cumulative.probs[1,] / cumulative.probs[1,6]
cumulative.probs[2,] <- cumulative.probs[2,] / cumulative.probs[2,6]
round(cumulative.probs, digits = 3)*100

#Percentage difference
round((cumulative.probs[2,] - cumulative.probs[1,]) / cumulative.probs[1,], digits = 2)*100

#Generate permutations for randomization inference
n.perms <- 100000
perms <- replicate(n.perms, genperm())

#Met people at each rank
for(i in 5:1){
	print(i)
	met.this.high <- as.numeric(data$staffrank >= i)
	print(ri(met.this.high, data$treat.donor, perms))
}

#Ordered probit
library(MASS)
staffrankfactor <- ordered(data$staffrank)
get.ll <- function(treat) logLik(polr(staffrankfactor ~ treat, method = "probit"))[1]
ll.actual <- get.ll(data$treat.donor)
ll.dist.sharp.null <- replicate(5000000, get.ll(genperm()))
mean(ll.actual <= ll.dist.sharp.null) #pvalue